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Superconducting circuits with Josephson junc- 
tions ' are promising candidates for developing 
future quantum technologies. Of particular in- 
terest is to use these circuits to study effects that 
typically occur in complex condensed-matter sys- 
tems. Here, we employ a superconducting quan- 
tum bit (qubit) — a transmon — to carry out an 
analog simulation of motional averaging, a phe- 
nomenon initially observed in nuclear magnetic 
resonance (NMR) spectroscopy ' . To realize this 
effect, the flux bias of the transmon is modulated 
by a controllable pseudo-random telegraph noise, 
resulting in stochastic jumping of the energy sep- 
aration between two discrete values. When the 
jumping is faster than a dynamical threshold set 
by the frequency displacement of the levels, the 
two separated spectral lines merge into a single 
narrow-width, motional-averaged line. With si- 
nusoidal modulation a complex pattern of addi- 
tional sidebands is observed. We demonstrate ex- 
perimentally that the modulated system remains 
quantum coherent, with modified transition fre- 
quencies, Rabi couplings, and dephasing rates. 
These results represent the first steps towards 
more advanced quantum simulations using arti- 
ficial atoms. 

The ability to resolve energy variations AE occurring 
in a time interval At is fundamentally limited by the 
energy-time uncertainty relation AEAt > h. Consider 
for example a system of spin- 1/2 particles filling a porous 
material (Fig. la) in an external magnetic field; then, 
if the pore wall is paramagnetic, the particles near it 
will experience a modified local magnetic field. As a re- 
sult, static particles will produce two spectral peaks in 
the NMR spectrum at frequencies ui and Ld2 (Fig- lb). 
In contrast, particles moving swiftly back and forth be- 
tween the two regions on timescales At shorter than 
h{AE)~^ = {uj2 — uJi)~^ are not able to discriminate 
between the two energy values. Then, as the particles 
move faster, the outcome in the spectroscopy is not sim- 
ply a continuous broadening and overlapping of the two 
peaks. Instead, a new peak emerges at the average fre- 
quency ljq (motional averaging) with a width smaller 
than the energy separation AE (motional narrowing). In 
atomic ensembles and condensed-matter systems, this oc- 
curs via fast variations of the chemical potential, molec- 
ular conformation, effective magnetic fields, lattice vi- 
brations, micro-electric fields producing ac-Stark shifts 
etc. Closely related phenomena are the Dyakonov- 



Perel effect , the Dicke line narrowing of Doppler spec- 
tra , and the quantum Zeno effect . 

In this paper we report the observation of motional 
averaging and narrowing in a single artificial atom (a 
qubit), with a simulated fast-fluctuating environment un- 
der direct experimental control. This should be con- 
trasted with the typical situation in condensed-matter 
systems, where one has a large number of particles and 
the experimentalist can only indirectly attempt to change 
the fluctuation rate, typically by modifying a thermody- 
namic function of state such as temperature or pressure. 

Our device is a circuit QED system' - , consisting of a 
tunable transmon coupled to a quarter-wavelength copla- 
nar waveguide resonator. The circuit schematic and the 
simplified experimental set-up are shown in Fig. Id. The 
applied bias flux $ through the qubit loop consists of a 
constant part $dc and a time-dependent part (t) with 
an amplitude much smaller than ^dc- The Hamiltonian 
of the transmon qubit ' is 

H{t) = ^[u;o+am^, (1) 

where ljq = i:jp[cos^(7r$dc/$o) + sin^(7r$dc/<&o)]' - 
Ec/h. Here, $o — h/2e is the magnetic flux quantum, 
the Josephson plasma frequency is Wp = •^/SEqEj /h, 
^x,y,z denote the Pauli matrices, the single-electron 
charging energy is Eq = e^/2(CG + Cb) « /i x 0.35 GHz, 
Ej = Eji + Ej2 ~ 24i?c is the maximum Josephson en- 
ergy, and d = {Eji — Ej2)/Ej « 0.11 denotes the junction 
asymmetry. We choose $dc so that wo/27r = 2.62 GHz, 
which is far detuned from the bare resonator frequency 
Wr/27r — 3.795 GHz, allowing the dispersive measurement 
of the qubit through the resonator by standard homodyne 
and heterodyne techniques . 

The time-dependent part h^{t) of the energy splitting, 
determined by <I'ac(0: is controlled by an arbitrary wave- 
form generator via a fast flux line. The random telegraph 
noise (RTN) is realized by feeding pseudo-random rect- 
angular pulses to the on-chip flux bias coil. Ideally, the 
dynamics of h£^{t) is a stationary, dichotomous Marko- 
vian process, characterized by an average jumping rate 
X and symmetrical dwellings at frequency values of ±^ 
(see Fig. Ic). The number of jumps is a Poisson process 
with the probability P„(i) = {xt)^e~^* /nl for exactly n 
jumps within a time interval t. This Poissonian process 
simulates the temporal variations causing motional aver- 
aging in atomic ensembles and condensed-matter systems 
(Fig. 1). 

In order to find the effective transition energies and 
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FIG. 1. Schematics of motional averaging and simplified experimental set-up. a, Illustration of a single pore filled 
with a medium consisting of spin- 1/2 particles that are either static or moving. White and grey areas denote different magnetic 
environments, producing the energy splittings ftcji and hijj2, respectively (see b,c). b, The schematic spectrum of the motional 
averaging. The solid lines denote the spectrum of static particles, while the dashed line corresponds to moving particles, c, 
A sample trajectory of the energy splitting h[Luo + in equation (1). d, Schematic diagram of the sample device and the 
simplified experimental set-up. Two asymmetric Josephson junctions (-Eji and Ej2) are capacitively shunted (Cb) and coupled 
(Cg) to a resonator. The time-dependent transition frequency of the qubit is controlled with fast fiux pulses on top of a dc 
bias (/flux = /dc + /pulse), generated by an arbitrary waveform generator. A vector network analyzer (VNA) is used for probing 
the resonator at ojprobc ~ ^^i- One microwave source (uj) is used for driving the qubit. Additional details are shown in Fig. SI 
of Supplementary. 



decoherence rates of the modulated system, we calcu- 
late the absorption spectrum'' '. For a qubit subjected 
to RTN fluctuations as in equation (1), we exploit the 
quantum regression theorem to find the spectrum 
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J{u,-u:o)r-r2\T\ /^-i/J-J(t)dt^ 



dr, (2) 



where the expectation value is taken over noise realiza- 
tions and T2 = + Ti/2 is the decoherence rate in the 
absence of modulation. For our sample r2 ~ 27r x 3 MHz 
(determined in an independent measurement). For the 
Poisson process £_{t), Eq. (2) gives (see details in 
Supplementary) 
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with r2 = r2 + 2x and uj = uj — ujq. In our measurement, 
we have recorded the population of the qubit 's excited 
state, as shown in Fig. 2. The qubit is excited with a 
transverse drive 



^{t) = hg cos{Lut)ax 



(4) 



The simulation of the occupation probability presented in 
Fig. 2a (see Methods and Supplementary) includes also 
the effect of power broadening due to strong driving am- 
plitudes (5 > ri,r2). 

Two intuitively-appealing limits result from Eq. (.3) 
(see also Fig. 2) at the opposite sides of a dynamical 
threshold defined by x = the case of slow jump- 

ing, X <C ^-iid resolvable energy variations, AEAt ^ 
f^/x ^ the qubit absorbs energy at ujq ± ^ with the 
total decoherence rate r2 = r2 +x- The correlation time 
of the displacement process h£^{t) is = (2x)~"'^ and the 
linewidth (full width at half maximum 2r2) broadens by 
the amount of the reduced mean life-time of a qubit ex- 
citation (exchange broadening'). In contrast, for fast 



jumping processes, x ^ d ^-^id when the variations are 
not resolvable, AEAt ~ /i^/x ^ ^, the qubit absorbs en- 
ergy only at the frequency ujq with r'2 = r2 +£_'^/2x- The 
increase in decoherence rate by C^/2x can be related to 
the excursions of the accumulated phase u{t) — ^(t) dt 
[cf. equation (2)] for single noise realizations, especially 
to the diffusion coefficient of the process^ ' . Surprisingly, 
the averaged spectral line is well-localized around the 
mean value, since C^/2x <C 2^. The effect of motional 
narrowing can be best seen in Fig. 2c as a raise in the 
height of the center peak because the additional deco- 
herence ^2/2x decreases with increasing jumping rate x- 
When X is comparable with ^, there is a cross-over re- 
gion where absorption is reduced and the peak broadens 
due to enhanced decoherence (see Fig. 2). This is impor- 
tant for the improvement of qubit dephasing times since 
it implies that a longitudinally coupled two-level system 
(TLS) fluctuator is most poisonous when its internal dy- 
namics occurs approximately at the same frequency as 
the coupling to the qubit. 

A counter-intuitive aspect of motional averaging is that 
the system is at any time in either one of the states ±h^ 
and spends no time in-between. Yet spectroscopically it 
is not seen in the states ±/i^; instead it has a clear sig- 
nature in the middle. Also, the spectrum (3) decays as 
I/lu'^ far in the wings, showing a non-Lorentzian charac- 
ter that originates from the non-exponential decay of the 
correlator {a-{T)a+{Q))^'' . This is because the frequent 
(T^-changes have a similar effect to a continuous measure- 
ment of the system, slowing its dynamics in analogy with 
the quantum Zeno effect ' . 

To further explore these effects, we have used sinu- 
soidal waves to modulate the qubit energy splitting: 
h[ujQ + Scos{^lt)]. In Fig. 3a, b, the experimental data 
and the corresponding numerical simulation are shown. 
Resolved sidebands appear at w = ± kfl [k — 
0,1,2,...), and arc amplitude-modulated with Bessel 
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FIG. 2. Motional averaging, a, Numerical simulation of the occupation probability plotted in the — x plane with 
parameters 5/27r = 43 MHz, g/2n — 8 MHz, Ti/2n — 1 MHz, and r2/27r = 3 MHz. b, Measured excited state occupation 
probability Pe under the RTN modulation, c, Horizontal cuts of the measured occupation probability in b, denoted with 
white lines, at three different jumping rates (solid lines). For clarity, the curves are displaced vertically by 0.3. The dashed 
lines indicate the corresponding results of the numerical simulation in a. See Fig. S2 of Supplementary for the zero-detuning 
{uj = LUo) occupation probability at different ^. 



functions Jk{S/fl), where S and denote the modulation 
amplitude and frequency, respectively. The Hamiltonian 
_ff +_ffdiivc [equations (1) and (4)] can be transformed to a 
non-uniformly rotating frame'^'^'"" (see Supplementary). 
The theoretical prediction for the steady-state occupa- 
tion probability on the resolved sidebands > g > T2) 
is 
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which is compared to experimental data in Fig. 3c. These 
population oscillations can be alternatively understood 
as a photon-assisted version of the standard Landau- 
Zener-Stiickelberg interference (see Fig. S5 of Sup- 
plementary). Also, when the qubit drive is tuned close 
to the resonator frequency u ^ LOr and when the driv- 
ing amplitude g is large, the system can be seen as a 
realization of a quantum simulator of the ultrastrong 
coupling regime^'^. In our setup, the simulated coupling 
rate is estimated to reach over 10 % of the effective res- 
onator frequency, comparable to earlier reports" ' where 
the ultrastrong coupling was obtained by sample design. 
The sinusoidal modulation allows a different perspective 
on motional averaging. The RTN noise comprises many 
different modulation frequencies ft (Fig. Ic). Each fre- 
quency creates different sidebands {k = 1,2,3,.. .), which 
overlap and average away. At large fi's only the central 
band k = survives because the only nonvanishing Bessel 
function of zero argument is Jo{0) — 1. 

Finally, we show that it is possible to use the mod- 
ulated system as a new, photon-dressed qubit. Indeed, 
we can drive Rabi oscillations on the central band and 
on the resolved sidebands (see Fig. 4 and Fig. S3 in 
Supplementary). For fast RTN modulation, x ^ ^) 
and under the strong drive in equation (4), we con- 
struct (see Supplementary) a master equation describ- 
ing the Rabi-oscillation with the detuned Rabi frequency 
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C2 



^ and with the total decoherence 
rate Fj = r2 -I- ^ jlx- In the case of sinusoidal modula- 
tion, the Rabi frequencies on the sidebands are obtained 

^sg'§=g\U5m. 

To conclude, we have demonstrated the first realiza- 
tion of motional averaging and narrowing in a single- 
particle system, i.e., a superconducting transmon whose 
energy splitting was modulated between two predeter- 
mined values while maintaining random dwell times in 
these two states. We have also shown that the system 
preserves quantum coherence under modulation, with re- 
duced power broadening on the central band. Our work 
could open up several research directions. For example, 
a key limitation in solid-state based quantum processing 
of information is the decoherence due to the fluctuating 
TLSs in the dielectric layers fabricated on-chip- We 
anticipate that the dephasing times of the existing su- 
perconducting qubits can be dramatically improved if one 
can accelerate the dynamics of the longitudinally-coupled 
TLSs. The quantum coherence of the resolved sidebands 
and of the central band suggest that low-frequency mod- 
ulation offers an additional tool for implementing quan- 
tum gates, in analogy with the use of vibrational side- 
bands in the ion trap computers- ' . Finally, the recent 
progress in the field of nanomechanical oscillators makes 
possible the study of frequency jumps in nanomechanical 
resonators , which is predicted to be accompanied by 
squeezing . Our work paves the way for further simula- 
tions of quantum coherence phenomena using supercon- 
ducting quantum circuits. 

Methods Numerical simulations resulting in the occu- 
pation probability of Figs. 2a, c, 3b, 4a, and S2 and S5c in 
Supplementary are done by solving the master equation 
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FIG. 3. Occupation probabilities of the excited state of the qubit under sinusoidal modulations, a, Measured 
excited state occupation probability with amplitude S/2n — 250 MHz. b, Simulation with Ti/2n = 1 MHz, T2/2n = 3 MHz, 
and g/2n = 20 MHz. c, The cut along the dashed line in a,b at Q/2n = 210 MHz. The (black) circles are the experimental 
data. The (blue) line denotes the analytical occupation probability obtained from equation (5) by taking k — —1,0, 1, which 
overlaps perfectly with the simulation in the parameter regime, Q ^ g > T2- For fe = (0-band), the linewidth is narrower 
than the bare linewidth without modulation (see Fig. S4 of Supplementary). 
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FIG. 4. Rabi oscillations with modulations, a, Rabi oscillations versus driving frequency with RTN modulation: measure- 
ment and numerical simulation with loo/2tt = 2.610 GHz, ^/2-k « 95 MHz, x ~ 600 MHz, g/2-K = 91 MHz, ri/27r = 1 MHz, 
and r2/27r = 3 MHz. The corresponding Rabi oscillations in the case of sinusoidal modulation are shown in Fig. S3 of Supple- 
mentary, b, The measured Rabi frequency of the random and sinusoidal modulation g and goff, respectively, as a function of 
Rabi-pulse amplitude (same scaling in the both panels). The solid and dashed lines are fits to the linear portion of the data, 
c, Rabi frequencies versus amplitude of the sinusoidal modulation. The horizontal axis is the dimensionless scaled modulation 
amplitude 5/Q, which is also the argument for the effective Rabi frequency g'^'^ — g\Jk{5/i^)\ plotted with g/2ir = 60 MHz 
and color coding for k: (red dotted) and 1 (blue dashed). The (red) circles and (blue) triangles indicate the measured Rabi 
frequencies at the centre of 0-band and 1-band, respectively. 



in the frame rotating at w and applying the rotating wave 
approximation (the pure dephasing rate = — ri/2). 
For the random modulation (Fig. 2a, c and 4a), we ap- 
ply the method of quantum trajectories^ ' and the imper- 
fections of the experimentally realized wave forms [the 
raising/falling time (~ 2.0 ns) and the sampling time 
(0.5 ns)] are included in the simulation. 
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I. MEASUREMENT SET-UP 

The electronic measurement set-up at room temperature is illustrated in Fig. SI. For qubit spectroscopy with 
RTN modulations, the dc flux bias and RTN modulations are generated by an Agilent 81150A arbitrary waveform 
generator (AWG). The qubit driving signal from an Agilent E8257D analogue signal generator and the CPW cavity 
probe signal from an Agilent N5230C PNA-L network analyzer are combined together by a Mini-Circuits ZFSC-2-10G 
power splitter /combiner and sent to the cavity input line of the dilution refrigerator. The signal from the cavity output 
line of the dilution refrigerator is amplified and detected by a PNA-L network analyzer. For sinusoidal modulations, 
the tones are generated by a R&S SMR27 microwave signal generator and they are added to the dc bias generated by 
the Agilent 81150A via a bias tee (not shown). Between the radio frequency instruments (analog signal generator and 
network analyzer) and their corresponding lines, dc blocks (DCBs) are used for breaking possible ground loops. All 
instruments are synchronized with a SRS FS725 Rubidium frequency standard (not shown in the figure). 

For Rabi oscillation measurements, high ON/OFF ratio Rabi pulses are generated by mixing a continuous microwave 
signal from the Agilent E8257D with rectangular pulses from a SRS DG645 digital delay generator via two identical 
Mini-Circuits frequency mixers. The Agilent N5230C PNA-L is used as a signal generator. Its output signal is split 
into two parts. One part is used for generating measurement pulses in a similar fashion to the Rabi pulses; the other 
part acts as a local oscillator (LO) signal to mix the output signal from the CPW cavity down via a Marki IQ mixer. 
The IQ data is filtered, amplified by a SRS SR445A pre-amplifier, and digitised by an Agilent U1082A-001 digitiser. 

Measurement controlling and data processing are done by MATLAB running on a measurement computer. The 
communication between the measurement computer and the instruments is realized through IEEE-488 GPIB buses. 
To generate RTN pulses, we use MATLAB's internal Poisson random number generator poissrnd(X) to obtain binary 
RTN sequences, and load the sequences into the Agilent 81150A AWG. Each binary RTN sequence consists of 50000 
data points and around 5000 random jumps in average. The mean jumping rate x (of AWG's output) is modulated by 
changing the clock frequency v of the AWG: x ~ 5000 x v. To verify this relation between x £md v, we observe the 
RTN sequences at different i^'s by a fast oscilloscope (10 GS/s), count the number n of edges (jumping events) for 
certain period of time t, and calculate the real mean jumping rate by its original definition, x = ^I't- As an example, a 
0.5 /iS long RTN sequence with estimated mean jumping rate x ~ 5000 x 10 KHz = 50 MHz is shown in Fig. Sic. 
25 jumps are counted during 0.5 /is, which gives x = 50 MHz. As long as < 100 KHz, the formula x ~ 5000 x u 
gives a good estimation of the mean jumping rate. For v > 100 KHz, the real mean jumping rate is smaller than the 
estimated one, due to the intrinsic ~ 2 ns rising/falling time of the AWG (pulses shorter than 4 ns have a certain 
chance to merge together). 



II. THE RANDOM MODULATION 

We consider the transmon qubit 

whose energy splitting has a fixed mean at u)q/2-k = 2.62 GHz and the time-dependent part /i^(t) is modulated in 
as random telegraph noise (see the main text). We do not observe transitions to the higher transmon levels in the 
experiments with the current driving parameters (Fig. 2 of the main text) thus it is sufficient to include only the two 
lowest transmon eigenstates in the analytical and numerical analysis. A sample path of the time-dependent energy 
splitting is shown in Fig. Ic of the main text and in Fig. Sic. The jumping between the two 'states' /i(a;o ± 5) occurs 
at random times, that is, /i5(<) is a stochastic process. We consider a special case of a stochastic process, namely 
a Poisson process, which is characterized by its mean jumping frequency x (#jumps/time) and by the probability 
-Pn(i) = (xi)"'e~'^*/n! for exactly n jumps within a time interval. Another way to characterize the stochastic process is 
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FIG. SI. Electronic measurement set-up. a, Detailed diagram of room temperature electronics for qubit spectroscopy, b, 
Detailed diagram of room temperature electronics for Rabi oscillations, c, An example of RTN modulation sequence generated 
by the AWG at a clock frequency u = 10 KHz, taken by an oscilloscope. 
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to construct the autocorrelation function R{t), 

R{r) ^ + r)m^ = e J2(-irPu{\r\) = ^^e->^l-l ^ = C^e'^^l-I, (SI) 

which defines the correlation time of the displacement process: = l/2x- As a notational convention, we will use the 
notation {^{t + T)^{t))^ to denote an average over the fluctuation process ^{t) without the risk of confusion with the 
usual ensemble average. 



A. The qubit spectrum 



The absorption spectrum of a qubit is defined as'^ 

SH = ^ r e'"-(a_(r)a+(0))^ dr. (S2) 

The correlation function (a- (t)(t+(0)){ is calculated by tracing over the environmental degrees of freedom and averaging 
over the random telegraph noise (RTN) modulations. The former, denoted as (a-{T)a+{0)), can be evaluated by using 
the quantum regression theorem ' , which states that if the equations of motions for the expectation values of a set of 
system operators are expressed as 

^(i.(t))=^G,,(i)(i,(t)), (S3) 
j 

then for the correlations the following relations hold 

d 



dt 



(i,(t + r)ifc(i)) =^Gy(T)(4(t + r)ifc(t)), (S4) 



where = Tr and p denotes the density matrix of the system. This means explicitly that if we know 

the equation of motion for ((7_(t)), then we know it for the correlation ((T_(r)(T_|_(0)) as well; only the initial value 
may differ. From the master equation 

^=-^ [H{t)J{t)]+'^{2a.pa+-a+a.p-pa+a.)-'^[a,,[a,,p]], (S5) 
we now first build the corresponding equation (S3) for ((T_(i)): 

^(a_(t)) = Tr = (^-i^o - m - Ti - ^) (^-(i))- (S6) 

Then, by the quantum regression theorem applied to equations (S3)-(S4), we can write the corresponding equation for 
the correlation function ((T_(t)(3"+(0)) and solve it, 

(a_(r)o-+(0)) =e-''^°^-^=l^'e-'-''o"fMdt^ (gj^ 
where we have set the initial value ((T_ (0)(t_|_(0)) = 1. The average over the fluctuations is denoted formally as 

(a_(T)a+(0))5 =e-'"°"-^=l"l {^e-'-Z'o" '^'^^ , (S8) 
which gives the spectrum in equation (S2) 

5(a;) = — /"e'^^-^^'-r^l^l/e-'/o^fWdA dr (2) 
27r 7_oo \ / 5 

= ^ / e'("-"°)--r^l^l (/(r))^ dr. (S9) 
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The next step is to calculate the average over the RTN fluctuations in /(t) = exp(— i JJ^ ^(t) dt) to get the explicit 
expression for the spectrum " . We start by forming a differential equation /(t) = — if(T)/(T), which, after a single 
iteration, assumes the following form 



df{r) 



Jo 



ar)mf{t)dt. (SIO) 



dr 

Now, by exploiting the symmetry of the jumping process = and the factorization rule' 

iahmh) ■ . . e(in))e = (e(il)e(i2))c(C(i3)C(t4) . . . atn))i, (Sll) 

where ti < t2 < ■ ■ ■ < in, we write the differential equation for the average of the function /(r) 



dr 



R{r-t){f{t))^dt. (S12) 



Noteworthily, the displacement process fi^it) does not obey Gaussian statistics in spite of the identical autocorrelations, 
implying that the usual Gaussian cumulant expansion cannot be applied . In other words, the above factorization 
rule in equation (Slf) does not hold for a Gaussian displacement process , whose corresponding correlation functions 
of order higher than two (n > 2) vanish. With an application of Laplace transformation to the integro-differential 
equation (Sf2), we notice that the Laplace transform F{s) of (/(t))^ obeys a simple algebraic rule 

F{s) = — i-, (S13) 
s + R(s) 

expressed with the Laplace transform R{s) of the autocorrelation function R{t) (see Ref. for an alternative derivation). 

Up to this point, the consideration has been general for any Markovian two-state noise process. The Fourier transform 
of time-symmetric functions can be calculated in an elegant way by exploiting the identity F{lij) = 7r~^Re[-F(s = — ii^)]- 
By using the explicit expression for the autocorrelation function R{t) in equation (Sf), the Laplace transformation 
becomes R{s) = ^^(s + 2x)~^- Thus, the Fourier transform of the function (/(t))^ is 



1 2y£2 



The explicit expression in time domain for the expectation value (/(t))^ is 

X 



(/(^)>C = e- 



^__sinh ( a/x^ - C^kl ) +cosh( Vx^-C^kl 



(S15) 



which is obtained with the inverse Laplace (Fourier) transformation. Here, one clearly observes the change of the 
dynamics in (/(r))^ as the dynamical threshold, defined by x = is crossed. The function (/(t))^ is oscillatory with 
jumping rates X < but it becomes a decaying function when X > 

The spectrum in equation (3) of the main paper is then given by the convolution of F(lu) in equation (S14) with a 
standard u; = ujq centered Lorentzian lineshape Suj„.r2{^) with full width at half maximum 2r2, 

5(c^) = ^ /"e'("-"°)--r^l-l/e-'/o^«*)dt\ (2) 

Q ( 'wr ' 1 2xg^ + r2[(r2+2x)^ + (a;-^o)^+e^] 

We derive the asymptotic limits of slow and fast jumping by considering the convolution form in equation (3). By 
reducing F{uj) in equation (S14) at the limit x ^ C around u) = ±^ and at the limit x ^ C around uj = 0, one gets that 

F{uj) = ^- A F{oj) = - ^ 2, (S16) 

' ^ 27r(a;T0' + x' ^^2+(^|^j' 

respectively. Combining these with the Lorentzian lineshape S^jg^r^i'-^') with the help of convolution in equation (3), 
one obtains the total decoherence rates in the limits of slow and fast jumping Fg = r2 + x and F2 = F2 + ^ jlx-, 
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FIG. S2. Zero-detuning {w = wa) occupation probability Pe versus jumping rate. The horizontal axis is the scaled 
jumping rate x/(C/27r)^. The markers indicate measured Pe at £,/2-k = 96 MHz (blue squares), i,/2TT = 71 MHz (green triangles) 
and (,/2-K = 43 MHz (red circles). The solid lines are the corresponding simulated Pe with the same driving and decoherence 
parameters as in Fig. 2a of the main text. 

respectively. Note that our convention for the units of frequency is that the jumping rate x is measured in frequency 
(Hz) and the other quantities, such as the jumping amplitude ^ or the decoherence rate in angular frequency 
(rad/s). 

In Fig. S2 we present the measured zero-detuning (cj = t^p) occupation probability at different measured jumping 
amplitudes ^ together with the corresponding numerical simulations (see Methods of the main text). The data sets 
show similar gradually increasing behaviour, and by plotting them as a function of the scaled jumping rate x/?^ thsy all 
fall onto same curve. This universal behaviour is explained by the absorption (without decoherence) in equation (S14) 
at the zero detuning F[ijj = 0) = 1x1 ■ The absorption increases linearly in xl^ until qubit saturation (reaching the 
steady state occupation probability = 0.5 by the moderately strong drive g > T2) for x/('?/27r)^ > 0.03 (MHz)~^. 



B. Master equation for a driven qubit 

When the qubit under RTN modulation of the energy splitting is driven with a relatively strong transverse (Rabi) 
drive, g > T2, equation (3) is not sufficient to describe the absorption since it does not include the saturation effects 
(power broadening) originating in the strong transverse driving. To overcome this limitation, we build up a master 
equation where the longitudinal RTN modulation is taken into account exactly in the qubit dynamics. We start by 
writing the total Hamiltonian (see equations (1) and (4) of the main paper) in the frame rotating with the driving 
frequency cj and assuming the rotating wave approximation (RWA), 

^tot = -ffq + H^{t) +Hr='^ [{UJO - ^)Oz + ff^^] + \i{i)^z + ^r, (S17) 

where ffq denotes the driven qubit, is the Hamiltonian of the reservoir (environment) coupling to qubit via H^{t). 
Now, we assume that the displacement process h^{t) is caused by coupling to the environment (bath). The intention is 
to trace the bath out and get equations of motion for the density matrix p of the qubit only. 

In the interaction picture (denoted with prime), the equation of motion for the density matrix w of the total system 
reads 

fw'{t) = ^[H'^(t),w'{t)]. (S18) 

By iterating (S18) and taking the average over fluctuations in H^{t) (trace over the bath Hj.), the density matrix p' 
becomes 

P=PiO)+p^[^yiy'^i£ dr2... dr„([i?^(n),[^^(r2),...[J?^(r„),p'(0)]]])^. (S19) 

Now, we adopt the notation L'^{t) = [H'^{t), (■)] for the commutator, which can be also understood as a 3 x 3 matrix if 
the density operator p' is seen as a 3 x 1 vector. By assuming the statistical properties described by the factorization 
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rule in equation (Sll), the above iteration can be truncated'^', such that we obtain the master equation 



P'(t) = r {^m'^ir)) p'(t) dr. (S20) 




dr " Jo \ 'U 

Transforming back to the Schrodinger picture results in 

^p(t) = ^L,p(t) - ^ ^* (Le(Oe-'^'(*--)L^(r))^ p(r) dr, (S21) 

where Lqp{t) = [Hq, p{t)] and e~'^«(*~^' is the time evolution operator, that is, the formal solution of the master 
equation p{t) = Lqp{t)/ih. Using now the vector representation for the density matrix p = {pz, Peg, Pge)'^ , where 
Pz = Pee ~ Pgg, wB have the following matrices corresponding to the Hamiltonian in equation (SI 7) 

/O 
\0 

The subsequent equation needs only information on the autocorrelation R{t) of the jumping process ^{t) in equation (SI) 
and the vector-matrix algebra, 

^ - rn'.P(^) - ^ / (L,(.)e-M-^)L,(.))^p(.)d.. (S22) 

To calculate the matrix exponent exp[— iLg(t — t)], we do a transformation to the eigenbasis of L^, where the matrix 
exponent has only entries 1 and exp[^i^ {ljq — w)^ + fl^(^ ~ t)] on its diagonal. Then we transform the master 
equation back to calculation basis, add the standard decoherence rates ^12, and express it using the Bloch pseudo-spin 
representation 

S^it) = [12 X S{t)l - ^ R{t) (^cos{giT)S,{t - r) + sin(3dr)5^(i - r) j dr - r2S',(0, (S23) 

Sy(t) = [n X S{t)], - f R{r) ( g^ + (^-'^o)^cos(gar) ^^^^ _ _ ,^{g^r)SS " r)) dr - r,5,(i), 

Jo \ 9d 9d J 

(S24) 

Sz{t) = [n X S(t)], - Vr{Sz(t) - 5,0), (S25) 



where 11 = (g, 0, — and = \J (cjp — tj)^ + denotes the detuned Rabi-frequency. 

In general, the correlation time r^ = l/2x of the displacement process ^(t) can have values comparable to the 
system dynamics (Rabi frequency g is of the order of 27r x 10 MHz), since in our measurements the jumping rate x 
varies between 50 MHz and 600 MHz. In general, the Markov approximation [short correlations in R(j)\ is not valid 
for solving equations (S23)-(S25). The Markov approximation can be used only, when the autocorrelation -R(r) has 
non-negligible values near r = 0. This is the fast-jumping case, x ^ in which situation the Bloch equations become 

SS) = X S(£)], - fra + 5,(t), (S26) 



2x 

Sy{t) = [12 X S(0]j, - (v^ + Sy{t), (S27) 

Sz{t) = [12 X S(i)]z - V^{Sz{t) - Szo), (S28) 

describing the Rabi-drive with the detuning a;o — a;, the driving amplitude g, the relaxation rate Fi, and the total 
decoherence rate FJ, = r2 + C^/2x- The time-dependent solution of equations (S26)-(S28) can be used to reproduce 
the measured Rabi oscillations (Fig. 4a of the main text). The numerical simulation (shown in Fig. 4a) includes the 
imperfections of the experimental wave forms (see Methods of the main text), giving very good quantitative agreement 
with the data. 
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FIG. S3. Rabi oscillations with sinusoidal modulation, a, Experimental data for the 0-band (centered at W()/2-k = 
2.62 GHz) and for the 1-band (centered at cjq + f2) taken with modulation amplitude (5/27r Ri 140 MHz and modulation frequency 
Q,/2tt = 120 MHz. b, Numerical simulation. The Rabi frequency without modulation is g/2-K = 60 MHz. The effective Rabi 
frequency (with detuning) for the fc-band is Peff""^' ~ \/ (t^o + kVl — lj)^ + g^J'^{S/Q,) in equation (S33). The decoherence 
parameters used for the numerical calculation are Ti/ln = I MHz and r2/27r = 3 MHz. 



III. SINUSOIDAL MODULATION 



We now consider a transmon qubit whose energy splitting is modulated simisoidally: h[ijjQ + Scoa{Qt)]. Even though 
n <C Wo, the qubit can be excited with the low- frequency signal in the presence of an additional high-frequency drive 
ui. This can be seen as a multi-photon transition process involving quanta from both fields [in equation (5) of the main 
text], or, as a photon-assisted Landau-Zener-Stiickelberg- interference^ (LZS) between the dressed states | J,, n) and 
|t,n-l). 



A. Multi-photon transition processes 



Let us first consider the multi-photon transition processes associated with the sinusoidal-modulated Hamiltonian, 
defined in equations (1) and (4) of the main text 



(S29) 



The qubit is driven both in the longitudinal (a^) and in the transverse (a^) direction. We show below how to eliminate 
the time-dependence from the longitudinal drive. After moving to a non-uniformly rotating frame with the unitary 
transformation ' ' ^ " 



U = exp 



- ( LUgt + — sin Qt ) a. 



(S30) 



the effective Hamiltonian is H' = U^'HU + ih{dtU^)U, and by using the Jacobi-Anger relations, we get 



hg 
2 



e"^"'a+ -Fh.c, 



(S31) 



By assuming that the transverse drive is close to a resonance, that is, a; w cjq ± Ml {k = 0, 1, 2 . . .) and that the 
resonances are resolvable fl > g > T2, we transform back with U = exp[i(tjJo + kQ — Lj)taz/'i\ and ignore all but the 
resonant terms, i.e., all the fast rotating terms (rotating wave approximation, RWA). The resulting RWA Hamiltonian 
reads 



h 



{uiQ + kQ,- u})az + gJk 



(S32) 



describing, in the Bloch spin representation, precessions around the vector il = {gJk (^) ,0,t<Jo + kO, — oj^ with the 
effective Rabi frequency 



(fc) 

9eS 



(S33) 
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FIG. S4. Excited-state occupation probabilities with and without modulation, a, The same as Fig. 4c in the main 
text, except for tlie red curve that indicates occupation probability without any modulation in equation (S34), with simulation 
parameters gjl-K = 20 MHz, ri/27r = 1 MHz, and r2/27r = 3 MHz. b, The blue curve denotes the measured occupation 
probability without any modulation and with the same driving amplitude as that in a; the red curve is the the same as the red 

curve in a. The linewidth of the 0-band (blue curve in a), -|- r2 [gJo (S/^)]'^ /Fi [see equation (-5)], is always smaller than 

the linewidth of the bare qubit (red curves in a, b), y/r'i + FagVri [see equation (S34)], since the Bessel function Jo{S/Q) < 1 
except at 5 = 0. 



Exactly at the multi-photon resonance u; = luq i: kfl (fc = 0, 1, 2 . . .), the effective Rabi frequency is \gJk{S/Q)\, and it 
is plotted in Fig. 4c of the main text. 

In the time domain, the measured Rabi oscillations are presented in Fig. S3 together with numerical simulations 
of Bloch equations exploiting the effective Rabi frequency in equation (S33). To find the steady state occupation 
probability in the presence of relaxation with rate Fi and decoherence with rate F2, we solve the Bloch equations 
analytically. In the RWA, that is when > p > F2 is satisfied, one is allowed to add up independent contributions 
from all the resonances (the resolved sidebands). The result for the steady state excited state occupation probability is 



■sp 2ri [g-^fc inJi 

,±-^ Fi + (^0 - + kQ)^ + ft [gJk m'' 



In the non-modulated case, the corresponding expression for the occupations is simply 

The occupation probabilities in equation (.5) and (S34) are compared with the experimental steady state occupation 
probabilities in Fig. S4. An interesting effect (Fig. S4b) is that the linewidth of the modulated qubit on the central band 
is always smaller than the linewidth on the qubit in the absence of the modulation due to reduced power broadening. 



B. Photon-assisted Landau-Zener-Stiickelberg interference 



The spectra seen in the experiment under sinusoidal modulation can be also interpreted as a photon-assisted 
Landau-Zener-Stirckelberg (LZS) effect. Note that in the absence of the driving field standard LZS processes are 
not possible: indeed the qubit energy separation is one order of magnitude higher than the modulation frequency, 
therefore the standard LZS probability is negligibly small. However, the system can still perform LZS-transitions by 
absorbing a photon from the driving field. This photon-assisted LZS-interference can be seen by transforming the 
Hamiltonian H + Hdiive from the equation (S29) into the frame rotating at ui around the a-axis [unitary transformation 
U = exp(— ia;(T2/2)], where it has exactly the same form (in the RWA) as that of LZS-interference " , namely 

H = ^[loq — uj + S cos{i}t)]a2 + ^ga^. (S35) 

As illustrated in Fig. S5a, LZS-transition events may occur when the modulation amplitude S is of the order of, or 
larger than, the detuning between the driving frequency and the qubit splitting lj — ujq. The LZS-transitions occur 
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FIG. S5. Photon-assisted Landau-Zener-Stiickelberg interference, a, Schematics of the photon-assisted LZS interfer- 
ence. The black diagonal lines denote the transversely dressed energy levels of the transmon |'|' / J,, Ji), shifted with photon 
quanta ±nhuj from the original energy levels of the transmon. The red curve shows the 27r/Q-periodic time-evolution of the 
energy splitting e = h[ujo + 5 cos(f2f)]. The blue curves are the adiabatic energy levels of Hamiltonian in equation (S35) that 
the system would follow in the absence of LZS-tunneling. b, Measured excited-state occupation probability as a function of 
the sinusoidal wave amplitude 5 (in dBm) added on top of dc flux bias, with fixed modulation frequency d/2n = 50 MHz. c, 
Simulation corresponding to panel b, with parameters ri/27r = 1 MHz, r2/27r = 3 MHz, and 5/27r = 20 MHz. 



between the transversely dressed states | J., n) and | 'l,n + 1), where n refers to the photon number of the transverse 
driving field. The phase difference of the two states gathered between the consecutive tunneling events leads to either 
constructive or destructive interference observed as maxima or minima in the occupation probability of the excited 
state, shown in Fig. S5b,c. The LZS-interference seen in our system corresponds to the so called fast passage limit '. In 
this limit, the expression for the occupation probability can be calculated analytically and the result agrees exactly 
with our equation (5). 

To prove experimentally that this picture is valid, we have scanned the modulation amplitude S and the driving 
frequency u) at fixed fl. We observe an interference pattern of the steady state occupation probability in good agreement 
with the theoretical prediction, see Fig. S5b,c. 
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